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Abstract 



We extend our previous first-principles theory for perovskite ferroelectric 
phase transitions to treat also antiferrodistortive phase transitions. Our ap- 
proach involves construction of a model Hamiltonian from a Taylor expansion, 
first-principles calculations to determine expansion parameters, and Monte 
Carlo simulations to study the resulting system. We apply this approach to 
three cubic perovskite compounds, SrTiOa, CaTiOs, and NaNbOs, that are 



O ' known to undergo antiferrodistortive phase transitions. We calculate their 

^ ■ transition sequences and transition temperatures at the experimental lattice 

constants. For SrTiOs, we find our results agree well with experiment. For 
^ ] more complicated compounds like CaTiOs and NaNbOa, which can have many 

OO I different structures with very similar energy, the agreement is somewhat less 

satisfactory. 
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I. INTRODUCTION 

Perovskite materials are of considerable interest both for fundamental reasons and for 
their many actual and potential technological applications. The great fascination of the cubic 
perovskite structure is that it can readily display a variety of structural phase transitions, 
ranging from non-polar antiferrodistortive (AFD) to ferroelectric (FE) and antiferroelectric 
(APE) in nature.El The competition between these different instabilities evidently plays itself 
out in a variety of ways, depending on the chemical species involved, leading to the unusual 
variety and richness of the observed structural phase diagrams. Por example, as temperature 
is reduced, BaTiOa undergoes a series of PE phase transitions, while SrTiOs has a single 
APD transition. More extreme examples are NaNbOs and BaZrOs; the former has a series 
of six transitions, while the latter stays cubic down to zero temperature. Another appealing 
property of these cubic perovskites is that all of the structural phase transitions involve only 
small distortions from the ideal cubic structure, the typical distortion being less than 5% 
of the lattice constant. This simplifies the theoretical treatment considerably. The ample 
experimental data on these compounds also provide many insights and opportunities for 
checking the accuracy of theoretical calculations. 

It is no wonder that there have been many theoretical attempts to study these 
compounds. Previous phenomenological model Hamiltonian approachesnia have largely 
been limited by oversimplification and ambiguities in interpretation of experiment, while 
empiricalcl and non-empirical pair-potential methodaZi have not offered high enough accuracy. 
Recently, advances in density-functional techniques have made possible first-principles inves- 
tigations of such perovskite compounds. Such calculations have proven capable of providing 
accurate structural properties and PE distortions for perovskites at zero temperature. ETEJ 

Recently, a thermodynamic theory based upon such first-principles calculations was de- 
veloped to study the finite-temperature properties of BaTi03,liiHi3 and predicted the cor- 
rect transition sequence and fairly accurate transition temperatures. This thermodynamic 
approach involves three steps: (i) constructing an effective Hamiltonian to describe the 
important degrees of freedom of the system; (ii) determining all the parameters of this effec- 
tive Hamiltonian from high-accuracy ab-initio LDA calculations; and (iii) performing Monte 
Carlo simulations to determine the phase transformation behavior of the resulting system. 
A similar approach was also successfully applied to PbTiOa by Rabe and Waghmare.ta 

The construction of the effective Hamiltonian is carried out in view of the special struc- 
ture properties of cubic perovskite compounds. At higher temperature, the cubic perovskite 
compounds ABO3 have a simple cubic structure with O atoms at the face centers and metal 
atoms A and B at the cube corner and body center, respectively. The two most common 
instabilities result from the softening of either a polar zone-center phonon mode, leading to a 
PE phase, or the softening of a non-polar zone-boundary mode involving rotations of oxygen 
octahedra, leading to an APD phase. (In some cases a zone-boundary polar mode may also 
occur, leading to an APE phase.) In our previous thermodynamic theory for BaTiOs, we 
assumed PE and strain distortions would be the only important degrees of freedom of the 
system. In other words, all other distortions are assumed to be much higher in energy. This 
is true for BaTiOs, but not true for cubic perovskites in general. As shown in our recent 
first-principle calculations,^ most cubic perovskite compounds may also undergo APD tran- 
sitions. To study these compounds, we need to extend our theory to include APD distortions 



among the low-energy distortions. This extended theory would also allow us to study the 
interaction between FE and AFD instabilities. However, because the AFD distortion is a 
zone-boundary distortion without a clear corresponding zone-center mode, the extension is 
not trivial. 

The rest of the paper is organized as follows. In Sec. 11, we go through the detailed pro- 
cedure for the construction of the effective Hamiltonian with the AFD distortion included. 
In Sec. Ill, we describe our first-principles calculations and the determination of the expan- 
sion parameters for the three compounds SrTiOa, CaTiOa, and NaNbOs. In Sec. IV, we 
report our calculated transition temperatures, phase sequences, and order parameters for 
those three compounds. We also identify the differences between the correlation functions 
of the FE and AFD local modes in SrTiOs. Sec. V concludes the paper. 

II. CONSTRUCTION OF THE HAMILTONIAN 

A. Local modes for AFD distortion 

In our previous development ,li3 we argued that the total energy of a cubic perovskite 
can be well approximated by a low-order Taylor expansion over all the relevant low-energy 
distortions, specifically FE distortions and strain. The FE distortions are represented by 
local modes, whose arrangement will reproduce the FE soft phonon modes throughout the 
Brillouin zone (BZ). To extend the theory to include the AFD distortions, we need to 
construct a new set of local modes to represent the lowest AFD modes over the whole BZ, 
or at least over the portion of the BZ where the energy change due to the AFD distortions 
is either negative, or positive but small. The AFD mode typically has the lowest energy at 
the zone-boundary R (0.5, 0.5, 0.5)27r/ao and M [(0.5, 0.5, 0)27r/ao, etc.] points, while near 
the zone center F the energy is very high. So it is necessary to choose local modes that will 
accurately reproduce the potentially soft modes in the vicinity of the R and M points. 

The rotation of an isolated oxygen octahedron can be represented by a pseudovector 
passing through its center. Assuming the origin of coordinates at the center of the octahe- 
dron, a pseudovector with polarization z involves displacements ±0.5aoy for oxygen atoms 
at (±ao/2,0,0) and displacements =F0.5aox for atoms at (0,±ao/2,0). Here oq is the lattice 
constant of the ideal cubic perovskite. In the case of the ABO3 perovskite crystal, we can 
represent octahedral rotation using pseudovectors sitting on the center of each octahedron, 
i.e., on the B atoms. However, neighboring octahedra share oxygen atoms, so that some 
continuity conditions would have to be imposed if we were to insist that the displacement 
of a given oxygen be consistently described by both neighboring pseudovectors. With such 
constraints the neighboring pseudovectors would no longer be independent of one another, 
leading to potential problems in the implementation of the Monte Carlo simulations. 

To avoid such problems, we simply construct a set of "virtual" pseudovectors aj which 
are independent of each other, and let the actual oxygen displacements be the superposition 
of the displacements that would result from these. To be precise, let a^ = a(Rj) denote the 
pseudovector centered on the B atom of unit cell i (position vector Rj), so that each oxygen 
atom is shared by two pseudovectors. The physical displacement of the oxygen atom shared 
by aj and slj is then given by 



Ar = yRij X (Hi - a^) , (1) 

where Hij = R, — Rj, Rjj = Rjj/|Rjj|. The AFD soft modes of interest at R and M 
are then easily represented by the corresponding pattern of pseudovectors. For example, 
choosing a(/x + my + nz) = (— )'"'"'"z reproduces one of the M-point modes polarized along 
z. (Here, x = aox, etc.) Other possible choices of the pattern of pseudovectors correspond 
to other modes which are probably higher in energy, but possibly still relevant for some 
materials. For example, choosing a(/x + ?7iy + nz) = (— )'y corresponds to an X-point mode 
that can be regarded as either of AFD character polarized along y, or of AFE character 
polarized along z. Finally, note that choosing a(Rj) constant (F-point arrangement) gives 
rise to no displacements whatever. 

In view of this last point, it is important realize that the aj themselves do not have 
direct physical meaning; only differences between adjacent a^'s are physical. Adding a con- 
stant to all aj's will not change the physical configuration of the system. So any physical 
distortions can be mapped to infinitely many pseudovector arrangements, but any pseu- 
dovector arrangement only corresponds to one specific physical distortion. Because only the 
pseudovector differences between sites have physical meaning, the Hamiltonian should be 
expanded in terms of these differences, not the pseudovectors themselves. Using this ap- 
proach, we reduce the number of degrees of freedom associated with oxygen displacements 
perpendicular to the 0-B bonds from six to three, and raise the symmetry of the system 
considerably. As a result, the Taylor expansion is significantly simplified. 

The two other low-energy distortions, the FE and elastic distortions, are treated as in Ref. 
|T2|. Briefly, for each unit cell of the ABO3 perovskite structure, we define a FE local-mode 
centered on the B site, and a displacement local mode centered on the A site. The former 
is chosen in such a way that a uniform superposition of FE local modes reproduces the soft 
TO mode obtained by diagonalizing the force-constant matrix at the Brillouin zone center. 
The quantities fj = f (Rj) and Uj are the vector amplitudes of the FE and translational 
local modes, respectively, in the ith unit cell. (Note the difference in notation between 



Ref. |T2| and this paper.Ej) Inhomogeneous strains r] are expressed in terms of differences 
between Uj in neighboring cells, and we add six extra degrees of freedom to describe the 
homogeneous strain. Thus, the total energy i?*°* depends on the set of variables {fj, a^, Uj} 
and the homogeneous strain, and is expanded in a Taylor series in terms of these quantities. 
The expansion terms can be divided into four kinds, those involving the FE local modes 
fj alone, the AFD modes a^ alone, the strain variable Uj alone, and the coupling between 
them, 

^tot ^ E'^iif}) + E^i{a}) + ^^({u}) + E'^\{i, a, u}) . (2) 

The part of the energy involving the FE local modes alone, E^{{f}), contains the on-site self 
energy, dipole-dipole interactions, and short-range residual interactions. Their forms have 
been given by equations in Ref. |T^: Eq. (3) in Sec. II. B, Eq. (7) in Sec. II. C, and Eq. (9) in 
Sec. II. D, respectively.Ej The energy due to Uj alone, E^{{u}), is just the elastic energy, and 
its form has been given by Eq. (11) of Ref. O. Also, the energy terms representing coupling 



between the fj and Uj have been given in Eq. (14) of Sec. II. F of Ref. [12. 



It remains to present here the energy terms involving solely the AFD local modes cij and 
those representing the coupling of the aj to the fj and Uj. Their expressions are presented 
in the following. 

B. AFD energy terms 

The AFD local modes a^ are nonpolar and involve no dipole moment, so long-range 
dipole-dipole interactions need not be considered, unlike for FE local modes. Recalling that 
only the differences of the a^ between neighboring sites are physical, it is not appropriate to 
separate energy contributions into on-site and inter-site interactions, as we did for the f(.ll3 
Instead, we separate the interaction into harmonic and higher-order contributions, 

E^{{8l}) = E^'^"™({a}) + ^^'""'^"^({a}) . (3) 

In principle, all the AFD energy terms should be expanded in terms of the Ar expressed 
through Eq. (|l|). For intersite interactions, this would become very complicated because of 
the low crystal symmetry at the O sites. However, for harmonic terms, the expression can 
be simplified by expansion in terms of the aj directly with certain conditions enforced. In 

this case, we can write 

E^'^-({a}) = 1^:1: 4,a/3«-%/5 ■ (4) 

ij af3 

Here, a and f3 denote Cartesian components, and Jlj^ap i^ a function of Rj^ and should decay 
very fast with increasing | Rj^ | . We need to impose conditions on the Jij^ap reflecting the fact 
that the dependence of the energy on the a^ is only through differences between neighboring 
sites [Eq. (|T])]. The appropriate conditions are 

E JLp = , (5) 

j G plane m 

where the sum is over sites j such that Rij^p = mao. This reflects the fact that if we make a 
change a^ —>■ a^ + cz for all the pseudovectors in an x-y plane at a distance of m unit cells 
away from the site i, the resulting the "force" on the pseudovector on site i should vanish. 
It can be shown that with these conditions enforced, the interaction energy is only related 
to pseudovector difference between adjacent sites. 

The description in terms of the aj directly makes it possible to simplify the interaction 
matrix J^^^p by symmetry, since a^ is centered on the high-symmetry sites. For a cubic 
lattice, we have 

on-site : J^,„^= 24'^a/3 , 
1st nn : 4,^= [jf + {j^ - jf )%,J5„^ , 
2nd nn : 4^^= [jf + {j^ - jf )%,J5„^ 

3rd nn : ^4^/3= Je^c^fs + j^dij,a0ij,i3il - S^p) , 

4th nn: 4,^^= J^e,,,J^p , (6) 



where %,q-=1 if Rij has a non-zero a component and otherwise. We include in-plane 
interactions {Ri^a = Rj,a) to 4th neighbor, since this kind of interaction is much stronger 
than other interactions. The conditions Eq. (^ can then be simplified to 

4 + 2j^ + 2j^ + 2j^ = , 

J^ + 4J3^ + 4J6^ = . (7) 

Thus, the complicated harmonic intersite interaction matrix for AFD local distortions can 
be determined from seven independent interaction parameters. 

The structural phase transition problem is intrinsically an anharmonic problem. Since 
the harmonic modes may be unstable, it is necessary to introduce higher order terms. For 
simplicity, we first only consider on-site anharmonic contributions associated with oxygen 
atoms. Because of the tetragonal symmetry on the O sites, the lowest anharmonic terms 
are of fourth order. Since each oxygen involves two nearest neighbor AFD pseudovectors, 
this quartic term will take the form 

^A,quart^|^|) =J2Y1 «^ {K(R0 - S(R^ + d)]' + [a.(RO - a,(R, + d)]^} 

i d=±x 
i d=±x 

-|- cyclic permutations . (8) 

Here, x = oox, and a^ and 7^ are parameters to be determined from first-principles calcu- 
lations. 

In our previous work on BaTiOs the intersite FE interactions have been expanded only 
up to harmonic order. For AFD interactions the corresponding approximation would be to 
truncate the interactions between the AFD-induced displacements of different oxygen atoms 
to harmonic order. [Such terms are already included in Eq. (^ .] We find this approximation 
to be satisfactory for those compounds with weak distortions, as in the case of BaTiOa or 
SrTiOa. For CaTiOs, the AFD distortion is very large and the transition temperature is 
around 2000K. In this case, we find it necessary to include more complicated anharmonic 
terms, such as third-order intersite interaction terms, for the AFD distortions. In fact, such 
terms turn out to be responsible for inducing a displacement component corresponding to 
an X-point phonon (with both O and Ca character) in CaTiOs. For NaNbOa, although the 
distortion is not as strong (the highest transition temperature is around 700K), there are 
many structures with very close free energy. We find that inclusion of the third-order AFD 
terms does have a noticeable effect for these compounds, so we include these third order 
interaction terms for CaTiOa and NaNbOs. 

We consider only those third-order interactions between AFD modes on two or three 
neighboring lattice sites. We can follow the treatment of the harmonic intersite interactions 
by listing all the possible interactions and using symmetry arguments to eliminate forbidden 
terms. Following this approach leads to three kind of terms, and we would need three more 
parameters to fully specify the Hamiltonian. Since the third-order terms are relatively weak 
and the exact determination of the three parameters is costly, we investigate the relations 
between these three parameters, and use a simple argument to combine the three terms to 
form a single new term with only one free parameter to determine. 

6 



The AFD interactions involve only the displacement of oxygen atoms. The strongest 
energy difference is associated with the distortion of oxygen octahedra, or the change of the 
length of a nearest-neighbor 0-0 bond A/. We can start by analyzing A/ for two oxygens 
at (0,0, a/2) and (a/2,0,0). We approximate the total-energy change as solely due to A/. 
Expanding it as a function of the rotation vectors a{i) up to the third order, we obtain the 
desired third-order intersite terms. Using the short-hand notations 

a' (a;) = a(Ri + x) - a(Ri) , 
a'(-a;) = a(Ri - x) - a(Ri) , (9) 

the cubic coupling term involving one 0-0 bond is 

Bs[a'y{z)-a'^{x)][a',{z)+a',{x)]'. (10) 

The other 11 nearest-neighbor 0-0 bonds will give rise to 11 other terms and the overall 
total-energy contribution can be expressed as E^'^^^ = B^J^iWi, where 

w. = {+Kiy)-<{zWy{z) + a'Mr 

-Ki-y) - a'x{z)][a'yiz) - a;(-t/)]2 

+Ki-y) - <i-m^yi-^) + <i-y)? 

+ cychc permutations } . (11) 

This assumption that A/ is solely responsible for the cubic intersite interactions significantly 
simplifies the energy expression and reduces the number of expansion parameters from three 
to one. 

C. Coupling energy 

There are three kind of coupling energy terms: those between FE and elasticity, between 
FE and AFD, and between AFD and elasticity. 

For simplicity, we consider only on-site couplings. The coupling between fj and Uj {E^~^) 
has been given by Eq. (14) in Sec. II. F of Ref. |12|.ll3 Here, we expressions for E^~^ and 

The coupling between elasticity and AFD modes at lowest order can be written as 

^^~^({a}, {r]}) = T^Yl 5/a/3x?7KR'i)[ai,a(x)ai,/3(x) + ai,a(-x)ai_/3(-x)] 

^ ilal3 

+ cyclic permutations , (13) 

where Oj ^(d) = aQ,(Rj + d) — aa(Rj) and r7;(Rj) is the six-component local strain tensor in 
Voigt notation {rji = en, rj^ = 2623). '7«(R«) can be expressed as a function of u, following 



Sec. II. F of Ref. |T^.E3 Biap^ is a high order coupling tensor. Because of the symmetry, there 
are only four independent coupling constants in Bia/3'y, 

^lyyx ^Izzx ^2xxy ^2zzy ^Zxxz -'^Syyz ; 

^2yyx J^3zzx -^Ixxy ^Zzzy ^Ixxz ^2yyz ; 

^3yyx ^2zzx ^3xxy ^Izzy ^2xxz ^lyyz ; 

J^iyzx ^izyx ^5xzy ^bzxy ^%xyz ^Qyxz ■ V / 

All other elements are zero. 

The lowest-order coupling between aj and fj is linear with both Hj and fj. It takes the 
form 

Ef-^ = E G,,yf^{Yii)[ay{Yl, + z) - ayiYii - z)] 

i 

+ cyclic permutations . (15) 

However, in the AFD state, this term is zero. To account correctly for the coupling between 
aj and fj, it is necessary to include higher-order terms. The lowest non-zero term in the 
AFD state is quadratic in both cij and fj. Defining Wj a, by 

Wi,x = o E [(^xi^i + d) - a^{Yii)] (16) 

° d=±y,±z 

and similarly for Wi^y and Wj^^, we can write 

^2 ~ Z^[ ^xxxxJi^x'^i^x + GxxyyJi^x\^i,y + '^i,z) 

i 

+ cyclic permutations ] . (17) 

(In principle a term Gxyxyfi,xfi,yWi^xWi^y could also be included, but for practical reasons we 
have not done so in this work.) In summary, up to the fourth-order terms, the coupling 
between aj and fj is expressed as 



III. FIRST-PRINCIPLES CALCULATIONS 

The expansion parameters in the model Hamiltonian can be obtained from a set of 
first-principles calculations. We use density-functional theory within the local density ap- 
proximation (LDA). The technical details and convergence tests of the calculations can be 
found in Ref. [1^. The use of Vanderbilt ultra-soft pseudopotentialst3 allows a low-energy 
plane-wave cutoff to be used for first-row elements, and also allows inclusion of semicore 
shells of the metal atoms. This makes high-accuracy large-scale calculations of materials 
involving oxygen and 3(i transition-metal atoms affordable. A generalized Kohn-Sham func- 
tional is directly minimized using a preconditioned conjugate-gradient method.li^Hi^y We 

8 



TABLE I. FE soft-mode eigenvectors for ABO^ cubic perovskites SrTiOs, CaTiOa, and 
NaNbOa- O || and O _L indicate oxygen displacement parallel and perpendicular to O — B bond, 
respectively. 

SrTiOg CaTiOs NaNbOg 

U 0.472 0.698 0.449 

(,B 0.612 0.330 0.625 

Coil -0.287 -0.157 -0.232 

Co± -0.400 -0.436 -0.421 



use a (6,6,6) Monkhorst-Pack k-point meshll^ for single-cell calculations (216 k-points in the 
full Brillouin zone), and the corresponding smaller sets of mapped k-points for supercell 
calculations. 

The calculation of expansion parameters related to the FE modes follows the procedure 
presented in Ref. |12|, Sec. III. The soft mode eigenvectors for SrTiOs, CaTiOs, and NaNbOs 
as calculated by King-Smith et ai, are summarized in Table |. The calculated expansion 
parameters for the FE modes are given in the top portion of Table ||. 

The calculation of the AFD expansion parameters follows a similar procedure as for the 
FE ones. The AFD eigenvector itself does not need to be calculated, since it is determined 
by symmetry. The LDA total energy vs. AFD distortion, with polarization along x and z, 
and at k-points X=(7r/a,0,0), M=(7r/a,7r/a,0), and R=(7r/a,7r/a,7r/a), are calculated. The 
arrangements of the AFD local modes are the same as for the FE-mode calculations as 
shown in Figs. 3(a)-(f) of Ref. |I^ However, the arrangements at the F point (Fig. 3(a) 



of Ref. |T2D and at the X point (Fig. 3(b) of Ref. |T2D involve no actual distortions. So the 
work reduces to four 10-atom cell calculations and two 20-atom cell calculations needed to 
determine ji, j^, j^, j^, Jq, and j^ + 2j^. The decomposition of j^ and j^ follows the 
same argument as for the FE case. It is difficult to perform sufficient LDA calculations to 
carry out the decomposition, and probably not very important to do so. Instead, we rely 
on the heuristic that the interaction between two AFD local modes should be minimal (in 
practice, zero) when they are so arranged that reversing one of them induces a minimum 
bond-length change. For the AFD case, this leads to Jq + 2jf = 0, allowing j^ and j^ to be 
obtained separately. The fourth-neighbor interaction parameter j^ is obtained from LDA 
calculations on a 15-atom supercell. 

As mentioned before, we need to include the effect of third-order intersite coupling in 
the effective Hamiltonian in some compounds having large AFD distortions. This kind of 
interaction generates a coupling among three distortions: an R-point mode with polarization 
(110), an M-point (7r/a,7r/a,0) mode polarized along (001), and an X-point (0,0,7r/a) mode 
polarized along (110). To determine the strength of this coupling, we carry out a calculation 
with the above R-point and M-point distortions frozen in (20-atom supercell), and calculate 
the forces of X-point character. We find that the force on the A-metal atom (Ca or Na) 
is non-zero, and opposite in sign to the force on the oxygen atom. This is in qualitative 
agreement with the experimentally observed displacements in the low-temperature phase, 
which are also opposite in sign. We then calculate the projection of these forces onto the 
X-point mode, under the simplifying assumption that the latter consists of equal-amplitude 



TABLE II. Expansion parameters of the Hamiltonian for SrTiOa, CaTiOs, and NaNbOs. En- 
ergies are in Hartree. FE focal-mode amplitudes are in units of lattice constant (a= 7.30a.u., 
7. 192a. u., and 7. 396a. u. for SrTiOa, CaTiOs, and NaNbOs, respectively); AFD local-mode ampli- 
tudes are in radians. 







SrTiOa 


CaTiOg 


NaNbOg 


FE on-site 


K2 


0.0559 


0.0240 


0.0679 




a 


0.150 


0.023 


0.168 




7 


-0.191 


-0.006 


-0.256 


FE intersite 


ii 


-0.02034 


-0.01186 


-0.02378 




J2 


0.04274 


0.02750 


0.03078 




is 


0.005722 


0.002040 


0.006460 




J4 


-0.003632 


-0.002886 


-0.005446 




J5 


0.004882 


0.001132 


0.004820 




J6 


0.001416 


0.000672 


0.002358 




J7 


0.000708 


0.000336 


0.001178 


FE dipole 


Z* 


8.783 


6.768 


9.179 




eoo 


5.24 


5.81 


4.96 


ADF harmonic 


1^2 


0.162238 


0.022244 


0.095852 




3t 


0.010526 


0.086972 


0.034884 




3^ 


0.000820 


0.000544 


0.002360 




3i 


-0.002782 


-0.001398 


-0.000272 




3t 


-0.105414 


-0.112050 


-0.097678 




3t 


0.009460 


0.010792 


0.009752 




3i 


0.002577 


0.001262 


-0.000318 




3^ 


-0.001288 


-0.000632 


0.000160 




3i 


0.013769 


0.013956 


0.014868 


AFD 3rd Order 


Bs 


— 


0.0056 


0.0029 


AFD 4th Order 


a^ 


0.05433 


0.04970 


0.03775 




^A 


0.04706 


0.02414 


0.04301 


Elastic 


Bii 


5.14 


5.15 


6.63 




Bl2 


1.38 


1.22 


0.96 




B44 


1.56 


1.29 


1.07 


FE-E coupling 


-Dlxx 


-1.41 


-0.59 


-1.71 




^lyy 


0.06 


0.06 


0.50 




^Ayz 


-0.11 


-0.10 


0.00 


AFD-E coupling 


^lyyx 


0.260 


0.234 


0.316 




■rA 

^2yyx 


-0.068 


-0.008 


0.026 




^Zyyx 


0.000 


-0.034 


0.031 




VI.A 
^Ayzx 


0.044 


0.040 


-0.041 


FE-AFD coupling 


Gxy 


0.0061 


-0.0001 


0.0014 




^xxxx 


0.53 


0.72 


0.35 




^xxyy 


0.11 


0.29 


0.06 
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out-of-phase motion of the two atoms. This projection determines our third-order intersite 
couphng parameter B^. 

The calculation of the anharmonic coefficient a^ is performed with an R-point distortion 
polarized along (001). Typically a set of eight calculations are performed for each compound. 
The resulting LDA total energy is fitted to a polynomial Eq + €20^ + 040^ using standard 
least-squares methods to extract k^ and a^. A single R-point calculation with polarization 
in the (111) direction is used to extract 7"^. 

The four parameters describing the coupling between AFD modes and elastic strains are 
obtained by performing four more 10-atom supercell calculations: at the M point, a || (001), 
with an isotopic strain rji = ri2 = rj^ = 0.01; at the M point, a || (001), with strain rji = 0.01; 
at the X point, a || (001), with strain 771 = 0.01; and at the R point, a || (111), with strain 
tja = rj^ = t]q = 0.01. Extra care has been taken to ensure cancellation of errors due to 
k-point sampling and basis-size differences for the different unit cells involved. 

Finally, the couplings between FE and AFD modes are determined as follows. The 
harmonic coupling G^y is determined by considering a geometry in which the primitive cell 
has been tripled along the x direction, and for which fx is non-zero in two primitive cells 
and tty is non-zero in the third. The anharmonic couplings Gxxxx and Gxxyy are determined 
from a series of calculations on a doubled-cell configuration in which f (F) || (100) and a(i?) 
II (100) or (010). 

The calculated parameters for all three compounds are listed in Table |I|. We note 
that the intersite interaction parameters between AFD local modes have a much stronger 
anisotropy than those between FE modes. For FE, the j's show no marked anisotropy. 
(Of course, when the Coulomb interaction is included, the actual interaction between FE 
modes are quite anisotropic). On the other hand, for the first-neighbor AFD couplings, j^ 
is more than one order of magnitude larger than j^, which is reasonable since j^ involves 
no distortion of oxygen octahedra. For second neighbors, jf is again much larger than the 
others, confirming that the distortion of oxygen octahedra dominates the energy for AFD 
distortions. This observation that the in-plane interaction parameters are much stronger 
than the out-of-plane ones is what prompted us to include also the fourth-neighbor in-plane 
AFD interactions in the effective Hamiltonian. 

IV. RESULTS 

After the expansion parameters have been determined from first-principles calculations, 
the finite-temperature properties of the compounds can be calculated using Metropolis 
Monte Carlo (MC) simulations.Ey The details of the MC simulations involving FE and 



elastic distortions have been described in Sec. IV of Ref. |T2|. With the AFD distortions 
included, the number of degrees of freedom is increased from 6 to 9 per unit cell. All the 
details are very similar, the main difference being that the AFD degrees of freedom introduce 
many more possibilities for modes which may go soft. The primary candidates for soft AFD 
modes are three modes at the R point and one mode at each of the three M points, all of 
them involving only rigid rotations of oxygen octahedra. We thus have to consider a rather 
complicated set of order parameters, and we anticipate that complex phases may form. 
The results for the three different compounds SrTiOs, CaTiOs, and NaNbOa will be 
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FIG. 1. (a) Intersite correlations of AFD local modes {az{0,0,0)az{x,y,z)) in SrTiOs at 
T=150K. The correlations in the x-z plane (y = 0, top) and the x-y plane (z = 0, bottom) are 
shown. Each small square represents one lattice site; the origin lies at the center. Grey scale denotes 
the degree of correlation: black for perfect correlation, white for zero correlation, (b) Same for 
FE local-mode correlations (/^(0, 0, 0)fz{x, y, z)). (c) Corresponding FE local-mode correlations in 
BaTiOa at r=320K. 

presented in the three following subsections. Because of the strong sensitivity of struc- 
tural phase transition temperatures to the lattice constant and the well-known ~ 1% LDA 
underestimate of lattice constants, we concentrate on presenting calculated transition tem- 
peratures and transition sequences at the experimental lattice constants. We thus implicitly 
apply a negative fictitious pressure to the simulation cell, as explained in Sec. IV of Ref. |T^. 

A. SrTiOs 



Thermodynamic properties for this compound have been calculated and published in 
Ref. |14|. A pressure Pq = — 5.4GPa is applied to restore the experimental lattice constant. 
A transition from the cubic phase to a tetragonal AFD structure at 130K, and two further 
FE phase transitions at 70K and lOK, were predicted. A later quantum path-integral MC 
simulation revealed that quantum fluctuations suppress the FE phases entirely, and reduce 
the AFD phase transition temperature to 110K.Eil This gives excellent agreement with ex- 
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periment, which reveals a single AFD phase transition at 105K,E3 and no unambiguous phase 
transition (but the presence of "quantum paraelectric" behavio]&c3) at low temperature. 
Our calculated pressure-temperature phase diagram showed that the FE and AFD insta- 
bilities have opposite trends with pressure, and FE and AFD instabilities tend to suppress 
each other. 

We have performed some further simulations to investigate the behavior of the AFD 
and FE local modes, and in particular the nature of the intersite correlations for FE and 
AFD local modes in the cubic phase but just above the phase transition temperature. The 
M-point AFD modes do not appear to be important for SrTiOs, so we focus on the two 
vector order parameters f (F) and a(R) associated with the zone-center FE modes and the 
zone-corner AFD modes. Since the order parameters are vectors, the correlation functions 
are second-rank tensors. These can be calculated from our simulations as 

Sa/six, y, z) = {vaixo, ^o, zo)vf3{xo + x,yo + y, Zq + z)) , (19) 

where the average is taken over all the sites xo,yo, zq in the MC simulation cell and over all 
MC sweeps t. Here the v^ denote the components of the FE or AFD order parameters [fa 
for FE, aQ,(— )'"'"™^"' for AFD, where Rj = /x -|- my + m\. We can get a good picture of 
the nature of the correlation by investigating the diagonal elements [a = (3) only. Since the 
three Cartesian directions are equivalent, it suffices to present Szz{x,y, z). 

In Fig. |l](a), we show the calculated AFD correlation function S^^^{x,y,z) in the two 
planes x-z and x-y for SrTiOa at T=150K, in the cubic phase but just above the AFD phase 
transition temperature of 130K. We find that the correlations are quite strong in x-y plane, 
with a correlation length of about three lattice constants. Along the z direction, even the 
nearest-neighbor vectors are almost completely uncorrelated. Thus, the shape of the "equal 
correlation surface" for AFD local modes is disc-like. This is easy to understand on the basis 
of the RUM picture.0 Since the AFD local modes involve a rotation of oxygen octahedra, 
and any distortion of the oxygen octahedra involves a large energy cost (as shown by the 
large magnitude of ji, J4, and jg), the AFD octahedral rotations about z correlate strongly 
in the x-y direction. On the other hand, the rigidness of the octahedra does not impose any 
relation between z-polarized AFD modes in different z planes (as reflected in the very small 
J2). Thus, the pancake-like correlation naturally results. 

Fig. |l](b) shows the corresponding FE correlation function S^^{x,y,z) for SrTiOs at 
T=150K (the FE phase transition occurs at 70K). Its behavior is just the reverse of the 
AFD modes, being strong along the z direction and weak in the x-y plane, and resulting 
in a needle-shaped "equal correlation surface." This behavior is a direct consequence of 
the anomalously large mode effective charges in the cubic perovskite compoundsE^ which 
strongly suppress the longitudinal FE fluctuations and leads to the strong correlation of fz 
in the z direction. On the other hand, the transverse FE modes can easily go soft, resulting 
in a short correlation length for /^ in the x-y plane. 

The above picture of the correlation functions for AFD and FE local modes are presum- 
ably quite general for the cubic perovskite compounds. For the FE modes, we decided to 
repeat the calculations for the case of BaTiOs, where the AFD instability does not intervene. 
We show in Fig. |l|(c) the correlation function calculated in the cubic phase at T=320K, about 
20K above the FE phase transition temperature. We can see the behavior of the correlations 
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FIG. 2. The probability distribution of the Cartesian component of the local-mode variable Ux 
for SrTiOs in the cubic phase at T=140K (dashed line) and T=150K (solid line). 

is again needle- like and even more pronounced than for SrTiOs, presumably because we are 
closer to the FE transition temperature. (An elongated supercell was used to accomodate 
the correlations in this case.) 

As was done for the FE modes in BaTi03,cJ it is revealing to compute the equilibrium 
distribution of one cartesian component of the AFD order parameter a(R) in the cubic phase 
just above the phase transition temperature in SrTiOa. This is shown in Fig. ^, where it 
can be seen that the distribution looks approximately Maxwellian. This is indicative of a 
transition having a character much closer to the displacive than to the order-disorder limit. 

B. CaTiOg 

CaTiOs is one of the more complicated perovskites. Experimentally, it is found to have 
two stable phases, an orthorhombic phase at lower T and a cubic phase above 1530K.c3 Some 
recent experiments suggest that the transition is to a highly disordered cubic phase.Ba The 
room-temperature orthorhombic phase has a very complicated structure with a 20-atom unit 
cell. The displacements of all the atoms away from their idea positions have been determined 
in Ref. ^. The refined structure as a function of temperature has also been determined 
recently using X-ray diffraction.El'EiJ This complicated structure can be decomposed into a 
simultaneous freezing in of three AFD modes: an R-point mode polarized along (110) with 
rotation angle 0.20 (angles in radians), an M-point mode polarized along (001) with rotation 
angle 0.14, and an X-point mode polarized along (110). The X-point mode involves not only 
the rotation of oxygen octahedra, but also an associated displacement of Ca atoms. The 
ratio of O and Ca displacement is about 1 : —3 and the oxygen octahedral rotation angle is 
only about 0.03. 

For such a complicated structure, even a complete first-principles determination of its 
T = structure would be very difficult. However, we can arrive at a partial understanding 
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FIG. 3. (a) Calculated AFD order parameters a{R) of CaTiOa vs. T at Pq 
Calculated FE and AFD order parameters of NaNbOa vs. T at Pq = — 4.3GPa. 
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of this structure as follows. Our calculations show that the reference cubic structure is 
unstable towards either the R-point or the M-point AFD mode individually (negative uj'^), 
whereas it is stable with respect to the X-point mode individually (positive u"^, no double- well 
behavior). In fact, the X-point mode involves a strong distortion of the oxygen octahedron, 
and thus is far from soft. However, the symmetry of the crystal is such that if both the 
R-point and M-point AFD mode distortions are already simultaneously present, then the Ca 
and O atoms experience forces in the pattern of the X-point mode, as a result of the cubic 
anharmonic interactions discussed in Sec. IIB. Thus, under these conditions the crystal would 
necessarily acquire some X-point mode distortion. We therefor conclude that the appearance 
of the X-point mode distortion must be the result of the third-order coupling between R-, 
M-, and X-point AFD modes. This is a major reason why we chose to include the third-order 
coupling in our model Hamiltonian. The strain coupling is expected to be important in the 
determination of the actual magnitude of the X-mode distortion, complicating the problem 
and making a complete quantitative LDA determination of the X-mode amplitude difficult. 
To obtain the phase transition sequence, we start the MC simulation at a high tem- 
perature (T > 2000K) and equilibrate the system for 10,000 MC sweeps. An isotropic 
pressure Pq = —ll.SGPa is imposed to restore the experimental lattice constant, and all 
the subsequent simulations are done under this pressure. The temperature is reduced in 
small steps (as small as lOK around the transition temperature) with 30,000 MC sweeps 
at each T to ensure equilibration. The order parameters are accumulated over the last 
20,000 MC sweeps, after checking that they do not vary significantly over this period. In 
our simulation, we find that except for the R-point AFD order parameters a(R), all other 
FE or AFD order parameters are zero throughout the simulation. So the only phase tran- 
sition we observe is associated with a(R). Fig. |l](a) shows the calculated order parameters 
a(i?) as a function of temperature. (What we actually plot are the averaged maximum, 
intermediate, and minimum absolute values of the order-parameter components.) At high 
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temperature (T > 1800K) the system is in the cubic structure, where all the components 
are zero. The material goes through a phase transition at 1750K, where all three order 
parameter components increase simultaneously. Our phase transition temperature is close 
to the experimentally measured 1530K, but we obtain the wrong low-temperature structure. 
Ours is rhombohedral with a 10-atom cell, instead of orthorhombic with a 20-atom cell as 
observed experimentally. 

The difference between our theoretical result and experiment is not quite as dramatic as 
it might seem. Both the structure and the structural energy are very similar for the R-point 
and the M-point AFD mode polarized along (001). Note that in the observed structure, the 
amplitude of the M-point mode is about v^/2 times the the amplitude of the R-point mode 
polarized in the (110) direction. Thus, we can say that the main difference between our 
calculated structure and the observed one is that one component of the R-point AFD mode 
is replaced by an M-point mode in the observed structure. As argued above, the additional 
presence of the X mode is just a result of third-order anharmonic coupling. In fact, we 
find that if we artificially increase the third-order coupling constant B^ by a factor of 5, we 
recover the experimental T = structure in our MC simulations. 

Clearly, the fact that the structure is strongly affected by relatively weak anharmonic 
intersite interactions makes the determination of the correct low-temperature phase very 
difficult in CaTiOa. It is possible that a more careful treatment of the cubic intersite inter- 
actions (for example, an independent determination of the coupling constants associated with 
all three of the cubic anharmonic invariants discussed in Sec. II. B, or three-site or further- 
neighbor terms) might bring a better agreement with experiment, although one should not 
rule out the possibility that neglect of quantum fluctuationg^ll or intrinsic limitations of the 
LDA might be at fault. 

C. NaNbOs 

Experimentally, NaNbOa is probably the most complex cubic perovskite known. The 
high-temperature phase is the simple prototype cubic structure as in the other cubic per- 
ovskites. Below 910K, a whole series of structural phase transitions has been found and at 
least six more phases have been identified. As the temperature decreases, the compound first 
goes through a cubic-tetragonal transition at 910K with freezing in of a(R) modes polarized 
along one axis. There are then three orthorhombic phases present in the temperature range 
845-638K, the most complicated having a unit cell containing 24 NaNbOs formula units. 
All of these phases can be regarded as given by rigid rotations of oxygen octahedra, accom- 
panied by small induced X-point distortions. From 638K down to at least 170K, NaNbOa 
is antiferroelectric with an orthorhombic unit cell containing eight formula units. At even 
lower temperature, the crystal has been reported to transform into either a rhombohedraS 
or monoclinico structure. 

The complexity of the structural phase-transition sequence suggests the presence of sev- 
eral competing structural instabilities with very similar free energies. In principle, all the 
distortions involved in the observed structures of NaNbOs are included in our model. How- 
ever, it is not realistic to expect that the calculated structural energies and free energies will 
be in exactly the right order, given the complexity of the problem and the level of accuracy 
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of current first-principles based approaches. Nevertheless, we believe a first-principles study 
of NaNbOa is still important in identifying the most prominent distortions, as well as for 
demonstrating the limitations of such approaches. 

The determination of the structure is done using MC simulations on a cubic 12 x 12 x 12 
simulation supercell. An isotropic pressure Pq = — 4.3GPa is imposed to restore the experi- 
mental lattice constant, and all the subsequent simulations are done under this pressure. We 
start the simulation at very high temperature and equilibrate. The temperature is reduced 
in small steps ranging from lOK to 50K depending on proximity to a phase transition. At 
each temperature step, at least 40,000 MC sweeps are used to ensure that equilibrium is 
reached. The order parameters are accumulated over the last 30,000 MC sweeps. 

The calculated averages of order parameters a(R) and f (F) are shown in Fig. [l|(b) as a 
function of temperature. All other modes are found to be zero throughout the simulation. 
As was the case for CaTiOa, the averaged maximum, intermediate, and minimum absolute- 
value components are plotted. At high temperature (T > 800K), the system is in the cubic 
structure with all the order-parameter components close to zero. As T decreases to about 
700K, one AFD component increases rapidly and becomes significantly non-zero, and the 
structure transforms from cubic to tetragonal. With further decrease of temperature, a 
second component became non-zero, indicating the occurrence of an orthorhombic phase. 
Below 560K, a third AFD component grows and the structure becomes rhombohedral. At 
very low temperature (below 50K), the simulation also apparently shows a sequence of three 
ferroelectric transitions, and the compound ends up in a rhombohedral ferroelectric structure 
at very low temperature. 

Our first cubic-tetragonal phase transition compares favorably with experiment; we ob- 
tain the correct structure and underestimate the transition temperature by only ~20%. In 
the orthorhombic phase, however, the calculated structure is much simpler than the ob- 
served one. Only one orthorhombic phase seems to occur in our simulation. However, 
Fig. [^(b) shows signs of fluctuations occurring in the vicinity of this phase (these fluctua- 
tions persist even if the number of MC sweeps is increased significantly). This indicates that 
the orthorhombic phase is not very stable, and may involve a mixing of different phases. 
Moreover, the transitions do not appear very distinct in Fig. |l](b) as a result of finite-size 
broadening, so increasing the lattice size may help to resolve the different phases. However, 
the computational load increases rapidly with increasing lattice size, and it becomes im- 
practical to carry out simulations at much larger size. Our inability to get the correct AFE 
phase at room temperature is probably the most significant failure of our approach. Our 
zero-temperature structure is ferroelectric, but since the FE phases occur only at such low 
temperatures, it is likely that quantum fluctuations would need to be included to determine 
the actual low-temperature structure .EiJ 

V. DISCUSSION AND CONCLUSION 

In this and previous studies, we perform a series of ab initio studies of the thermody- 
namic properties of perovskite compounds. Without introducing any adjustable parameters, 
we have calculated structural transition sequences, transition temperatures, phase diagrams, 
and other thermodynamic properties based on first-principles calculations. For compounds 



17 



with simple phase transitions, hke BaTiOa and SrTiOa, our calculated thermodynamic prop- 
erties agree very well with experiment observations. For more complicated compounds like 
CaTiOs and NaNbOs, our results are less satisfying. 

There are two major sources of errors, the inaccuracy of LDA calculations and the im- 
perfection of our models. Our LDA calculations have been carefully performed to avoid 
possible errors, and convergence has been carefully tested. As for the intrinsic accuracy of 
LDA, our calculated structural parameters and energies are within a few percent of experi- 
mental values. Although this is the usual high accuracy observed generally for the LDA, it 
is unfortunately not enough for truly accurate determination of the thermodynamic prop- 
erties of perovskites. For example, it is embarrassing that we are forced to choose between 
carrying out the calculations at the theoretical equilibrium lattice constant or the experi- 
mental lattice constant (negative fictitious pressure); this choice can affect phase transition 
temperature by ~100%. We regard this as being the most important probable source of 
error in our calculations. 

It is also possible to improve our model Hamiltonian. For example, our restricted as- 
sumption for the form of the third-order intersite interactions may be lifted, resulting in a 
significantly more complicated model Hamiltonian. Also, other higher-order terms can be 
included in the Hamiltonian. It would also be possible to include more degrees of freedom 
per cell, or include eigenvector information from more k-points of the Brillouin zone when 
defining the local-mode vectors, to treat other phonon excitations more accurately.^ How- 
ever, in view of the current accuracy of first-principles calculations, we are not sure that 
these modifications would dramatically improve our results. 

Finally, we emphasize that all of the MC simulations reported here treat the atomic mo- 
tion purely classically. As mentioned above, we have recently reported results of quantum 
path-integral MC simulations showing that quantum fluctuations of the atomic coordinates 
(i.e., zero-point motion) can shift transition temperatures by tens of degrees, and in some 
cases even eliminate delicate phases.^ Certainly this remains an important avenue of in- 
vestigation for CaTiOs and NaNbOs, but we nevertheless think it unlikely that inclusion 
of quantum fluctuations would immediately resolve the discrepancies with the experimental 
phase diagrams for these compounds. 

In conclusion, we have extended our previous first-principles theory for perovskite ferro- 
electric phase transitions to treat also antiferrodistortive transitions. We apply this approach 
to the three cubic perovskite compounds SrTiOs, CaTiOs, and NaNbOs, and calculate their 
thermodynamic properties including phase transition sequences and transition temperatures. 
For SrTiOs, our calculated results are in good agreement with experiment. For CaTiOs and 
NaNbOs, our calculated structural transitions have the correct general trend and the tran- 
sition temperatures are in rough agreement with experiment, but the calculated transition 
sequences are not correct in detail. We attribute this to the larger distortions and many 
multiple competing instabilities in these compounds. For SrTiOs, we also analyzed the in- 
tersite correlations for both FE and AFD local modes, finding needle-like and pancake-like 
correlations respectively for FE and AFD modes as expected on physical grounds. 
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